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OO Abstract 

o ■ 

. This paper establishes a relation between scale-free networks and Markov chains, 

' and proposes a computation framework for degree distributions of scale-free networks. 

. We first find that, under the BA model, the degree evolution of individual nodes in a 

' scale- free network follows some non- homogeneous Markov chains. Exploring the spe- 

. cial structure of these Markov chains, we are able to develop an efficient algorithm to 

Oh' compute the degree distribution numerically. The complexity of our algorithm is 0(t^), 

(-H , where t is the number of time steps for adding new nodes. Wc use three examples to 

demonstrate the computation procedure and compare the results with those from the 
existing methods. 
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1 Introduction 

Complex networks describe a wide range of practical systems of high technological, 
biological, and social importance [1,2]. For example, the Internet, the World Wide Web 
(WWW), biological cells and communities of scientists can all be described as complex 
networks. 

Erdos and Renyi [3] started the early studies of complex networks as random graphs 
in 1960. Many years later. Watts and Strogats [4]'s construction of the small- world 
network in 1998 represents an interesting development for the study of complex networks 
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in that it was motivated by observations of real system behaviors (e.g., Milgram's six- 
degree connectivity [5] ) . A common feature of the random graph and small- world models 
is that the degree distribution (the probability of finding a node with k connections) 
decays exponentially with the number of connections. However, empirical evidences 
from the Internet and WWW, among other complex networks, show a fundamentally 
different picture, i.e., the tail of the degree distribution follows a power law. This led 
to the introduction in 1999 of scalc-frcc networks by Albert, Barabaasi, and Jcong in 
their pioneering works [6-8], and the start of a new phase in the study of complex 
networks. Recent studies [9-20] are characterized by empirical observations of scale-free 
behaviors in various practical systems and investigations of the formation mechanisms of 
scale-free network. A number of important properties in scale-free networks have been 
identified, such as the small-world character, the emergence of hubs, and robustness 
and frangibility. These properties show that scale-free networks can play an important 
role in the understanding of many complex and important systems. 

Two general features can be observed in many real- world networks: successive addi- 
tions of new nodes and certain preference in linking to existing nodes. Albert, Barabaasi, 
and Jeong proposed two mechanisms to characterize the evolution of a scale-free net- 
work [7, 8]: the growth mechanism, starting from mo nodes, the network grows at a 
constant speed, i.e., adding one node at each time step and connecting to ?77,(m < mo) 
existing nodes; the preferential attachment, the chance that an existing node receives 
a connection from a new node is proportional to the number of connections it already 
has. The authors show that, under these two mechanisms, a network evolves into a 
stationary scale-free state. Its degree distribution follows a power law with the degree 
exponent 7 = 2.9 ± 0.1 from simulation analysis and 7 = 8 from the analytical result. 
These results are significant for complex networks and the two mechanisms become the 
first model, referred to as the BA model, by which large networks can self-organize 
into a stationary scale-free state. Empirical evidences show that in many networks, the 
number of edges grows faster than the number of nodes. This leads to the investigations 
of m- varying BA models, such as Dorogovtsev and Mendes [20]. 

Our research is mainly motivated by the following observation. While analytical 
solutions of the degree distribution for some simple models, such as the BA model, can 
usually be obtained, one has to resort to simulation for the degree distribution when the 
mechanisms in model become more complex. This may inhibit the further development 
of the theory on complex networks. In this paper, we propose an alternative approach. 
We first find that the degree evolution of a complex network can be characterized by a 
sequence of Markov chains. By carefully analyzing the structure of these Markov chains, 
we then develop an efficient numerical method to compute the degree distribution of 
complex network models. To show the feasibility and efficiency of our numerical method, 
we compute the degree distribution of the basic BA model and two of its variants. 

We organize the paper as follows. In the next section, we review some of the exist- 
ing methods for network degree distributions. We then use Markov chains to capture 
network dynamics. Exploring the special structure of the transition matrices of the 
Markov chains, we develop an efficient algorithm to compute the degree distribution 
asymptotically. We use this algorithm to compute the exponent of the degree distri- 
bution of the BA model. In Section 3, we compute the degree distributions of two 
m-varying BA models. We verify our approach by showing that our numerical results 
for the BA model and its variants match very closely to the existing results from the 
analytical and simulation approaches. We conclude the paper in Section 4 by pointing 
out some future research opportunities. 

2 A Markov chain-based numerical method 

With the preferential attachment mechanism of the BA model, the probability that 
node i receives a connection from an up coming new node is proportional to its own 
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degree fcj [7], i.e., 

m) = (1) 

Assuming continuity of ki{t) and treating n(A;i) as its rate of growth, ki{t) then satisfies 
the following dynamic equation [7, 8] 

Under the initial condition that ki{ti) = m, the solution of this equation leads to 

k,it) = m{^r , P=l (3) 

where U is the time when node i joins the network, and the degree distribution 

P{k) ~ 2m^fc-T , 7 = 3. (4) 

Here, f3 is called the dynamic exponent while 7 the degree exponent. 

The above simple analytical method is often rcfcrccd to as the continuum (mean 
field) theory. Similar power law results for the degree distribution are also obtained 
using different analytical methods by other authors. For example, with the master- 
equation approach [14], Dorogovtscv, Mcndcs and Samukhin treat the degree ki{t) of 
a node « at a fixed time t as a random variable. Thus its probability P{k,ti,t) for the 
BA model has the following relation: 

P(fc, U, t + l) = '^P{k - 1, ti, t) + {l- ^)P(fc, ti, t). (5) 

Let 

P{k,t)= ^'^ Y''' . (6) 

Assuming that the limit P{k) = limt_>oo -P(fc, i) exists and limt_>oo t + 1) — 

P{k,t)\ = (this is an additional condition), the degree distribution satisfies equation: 

2P{k) - 2Sk,m = {k- l)P{k - 1) - kP{k), (7) 

and the network degree distribution can be obtained as 

^^"'^ k{k + i){k + 2y 

Krapivsky, Redner and Leyvraz's rate-equation approach [15] focuses on the number 
Nk{t) of nodes with k edges at time t. For the BA model, Nk{t) is shown to satisfy 

dNk {k - l)Nk-i{t) - kNkjt) ^ - 

— e;^w — ^ 

Asymptotically, Nk{t) = tP{k) and kNk{t) = 2m,t, leading to equation (7). 

While the above methods handle simple models, such as the BA model, well, they 
do not, so far from the best of our knowledge, render analytical solutions for more 
complicated models. In this case, one can usually use simulation. While simulation is 
widely applicable, it is usually quite time consuming and may not be flexible enough 
for in-depth analysis of network behaviors. Here, we propose a different approach to 
capture the network dynamics. 

Consider the degree Ki{t) of node i at time t. Following the increase of t, the 
sequence {Ki{t),t = i,i + 1, ...} is, based on the preferential attachment mechanism of 
the BA model, a stochastic process with the state space fl = {171,171+ 1, ...}. Here and 
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below, wc use the upper case K to emphasize the fact that the degree sequence is a 
stochastic process. The attachment mechanism also indicates that the future evolution 
of the process is independent of the past history, given its current state; but it is time- 
dependent. This shows that the process {Ki{t)} is in fact a non-homogeneous Markov 
chain [21], with time-dependent transition probability 



Pkj{t + 1) = P{Kiit + 1) = J I K,{t) = k} 



1_ A 

2t' 

k_ 
2f 

0, 



j = k 
j = k + l 
otherwise 



(10) 



for k 



, m • 



and 



1, 



J 



0, j^k 



(11) 



fov k > m + 1 — i. Thus, the dynamics of a node from the time it joins the network is 
described by a non-homogeneous Markov chain and the whole network (excluding the 
original nodes) is completely described by t non-homogeneous Markov chains, where t is 
the time of the observation. Let Pi{t + 1) be the one-step transition probability matrix 
of node i at time t. We have, fov t = i,i + 1, ... 



m+i) 



2t 



2t 



1 m-t-1 m+1 

^ 2t 2t 



m+t- 
2t 



m+t- 
2t 



(12) 



Let fi{t) be the probability vector (distribution) of Ki{t) for a given t, and 



iS,T) 



t+l 



T 

E 

i=S 



Mt+1), P{k,t + l) 



F,f^\k-m + l) 
t-S + 1 



(13) 



Here, S and T are two fixed integers between 1 and t. Their meanings will be clear 
in the computation procedure later. The desired degree distribution of the network is 
then P{k) = limt^oo P{k, t + 1). 

Let us examine (13) to sec what is involved in computing the network degree distri- 
bution. It is clear that P{Ki{i) = k} = 1 if k = m and otherwise. We then have the 
initial probability vector 

/i(i) = (l,0,0,...) = ei (14) 

for any i. By density evolution of Markov chain, the f-|-l-step probability vector fi{t+l) 
is given by 



fi{t + l) = ei ■ Pi{i + 1) ■ Pi{i + 2) ■ ■ ■ Pi{t + 1) , t = i,i + l, 



(15) 



where the dots represent matrix multiplications. This, together with (13), shows that 

computing the degree distribution requires the multiplications and summations of an 
infinite number of infinite matrices. It is not realistic to expected any meaningful 
analytical solution from these computations. Even numerical computation seems un- 
manageable. Fortunately, our past experience in infinite matrix computations [22] with 
a rectangle-iterative algorithm guides us to explore the special structure of the one-step 
transition matrices. This leads to dramatically simplified matrix manipulations and a 
highly eSicient algorithm. 
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We note that while the transition matrices of consecutive nodes are different, their 
structural similarities lead to the following relations 

eiP,(i) = eiPi(t), z = 2,3,...; t = i + l,i + 2, ... (16) 

eiP^{t)P,{t + l) = eiPi{t)Pi{t + l), i = 2,3,...;t^i + l,i + 2,... (17) 
and in general 

eiP^{t)P,[t + l)---P,{t + s) = eiPi(i)Pi(< + 1) • • • Pi(i + s), (18) 
for i = 2, 3, t = i + 1, z + 2, ... and s = 2, 3, ... . Substituting the above relations into 

T T 

Ftli ^ = E /^(* + 1) = E + 1) • ^'(^ + 2) • ■ • ^^(* + 1)' (19) 

we obtain the following key relation 
F^+P - ((• ■ ieiPsiS + 1) + ei)P5(^ + 2) + • • •) + ei)Ps(T + 1) • • ■ Ps{t + 1). (20) 

The computation of Pj+i becomes very easy with (20). We start from the inner 
most bracket. After one multiplication and one summation, we obtain a row vector 
whose first two elements are nonzero. The second round of multiplication and sum- 
mation lead to a row vector with the first three elements being nonzero, and so on so 
forth. The final result is a row vector with the first (t — 5 + 1) elements being nonzero. 
An efficient algorithm can be developed to implement this procedure. Obviously, the 
complexity of the algorithm is O(t^). 

We plot the log — log curves for P{k, t) for some different m and t as shown in Figure 
1, and use the least square method to fit the exponent 7 and the coefhcient c of the 
power-law under the BA model. Table 1 lists the numerical results for different m 
and t values. We observe that the degree exponent is independent of m and the value 
matches those of simulation and the analytical solution with the mean field method. 
The coefficient of degree distribution c is between 2m^ and 2m{m+ 1), again matching 
the theoretical value from the mean field method. Furthermore, results for m = 3 show 
that the coefficient c is independent of t, i.e., the network is stationary. 




Figure 1: The degree distribution of the BA model 



In Figure n the three lines from left to right correspond to three cases: (1) m = 
1, t = 150,000; (2) m = 3, t ^ 100,000, 150,000, 200,000; (3) m = 5, t = 150,000. 
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The line in case (2) is the overlap of three hnes corresponding to three different t values. 
This shows that the distribution is stationary. The three lines of the three cases are 
parallel, which further shows the degree exponent of the BA model is independent of 
m. 

Table 1: Degree exponent and coefficients of the BA model 



parameter m 


time t 


exponent 7 


coefficient c 


1 


150000 


2.960830 


3.147515 


3 


100000 


2.989636 


21.79266 


3 


150000 


2.990032 


21.89667 


3 


200000 


2.980275 


21.01711 


5 


150000 


2.978891 


52.58130 



3 The degree distributions of m-varying BA models 

Our numerical approach is feasible and can be efficiently applied to more complex 
models. Since the number of edges grows faster than the number of nodes in many 
networks as shown by empirical evidences, we compute the degree distributions of two 
cases of the BA model with m-varying functions in this section. 

3.1 Power function 

Let the number of new links added in time step t be mt^, < 6 < 1, i.e., the new 
node t will link itself to m,t^ different nodes already present in the system. 

We note that after t time steps, this case leads to a random network with N = t+mo 
nodes and approximately Jg mx^dx links. Then, the total degree number of the system 
at time t is 

Y,k,^2j\x'dx = f^/^\ (21) 

Assuming continuity of ki{t), it then satisfies the following dynamic equation 

— =mt n(fc,) =mt^ = (22) 

Under the initial condition is ki{ti) = mtl, where U is the time when node i joins the 
network, we solve this equation and obtain 

k,{t)=mt<i{'-)'^=mt\'-f, /3=^. (23) 

Hence the degree distribution at time t 

P{k,t)^^mT-^t^k-^, (24) 

Here, z is called the non-stationary exponent. We note that this type of m-varying 
function was first discussed in [20] . 

We now construct the Markov chain for the degree sequence {Ki{t), t = i,i + 1, ...}. 
The state space is = {m^, mj -|- 1, ...}, where mj = m[i^]. At time t, the probability 
that an existing node i will connect with the new node is given by 
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Hence, the one-step transition probabilities are 
Pkjit + 1) = P{Ki{t + 1) = j I Ki{t) = k} -- 
for k = rrii, nii+t — i, and 



for k > TTii + t — i. The transition probabihty matrix is 




j = k 

j = k + l 
otherwise 



(26) 



(27) 



1 



rrtije + l) mi(0+l) 
2t 2t 

_ (mi + l)(e+l) (mi + l)(6>+l) 
^ 2t 2t 



(mi+t-i)(6i+l) (mi+t-»)(6>+l) 



2t 2t 
1 

(28) 

for f = z, i + 1, ... . 

We now provide the computation results when 9 = 0.2. We note that the structure 
of the transition matrices here is similar to that of (12). The difference is that now rrii 
is not a constant, in general, but a step function of i, as shown in Table 2. 



Table 2: Intervals of m,i keep constant 



time i 


32 


243 


1024 


3125 


7776 


16807 


32768 


59049 


100000 


161051 




2 


3 


4 


5 


6 


7 


8 


9 


10 


11 



Therefore, relations (16), (17) and (18) hold for each interval, e.g., the interval 
(243, 1023). Thus we obtain the following important result 

T-,{32,t) i;,(32,242) , ^,(243,1032) , , j;,(59049,99999) , ^,(100000,*) 

= + H — T ^t+1 + ■ v^y) 

Similarly, the initial probability distribution is fi{i) = (1,0,0, ...) = ei for any i. Thus 
the same algorithm based on (20) can be used to compute the degree distribution P{k, t) 
for this network. 

Prom the computation results, we plot the log — log curves for P{k,f) for some 
different m and t as shown in Figure 2. Wc also list some numerical results in the 
Table 3. From the figure and the table, it is clear that this network self-organizes into 
a non-stationary scale-free network, with the degree exponent 7 « 3.5. 

3.2 Logarithmic function 

Let the number of new links in time step t be m In t. 

We note that after t time steps, the model leads to a random network with N = t+mo 
nodes and approximately Jg m In xdx links. Then, the total degree number of the system 
at time t is ^ 

Y^kj !v2 j mlnxdx = 2mt{\nt - 1). (30) 
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Table 3: Numerical results of the power function case 



paramcler m 


time / 


(;xi)ou('Ul 


cocllicic'ut c 


1 


150000 


3.502938 


891.641 


3 


100000 


3.499978 


8213.46 


3 


150000 


3.502746 


10920.8 


3 


200000 


3.496971 


12300.2 


5 


150000 


3.503176 


37303.5 



The average degree of the system is A; = 2m{lnt — 1), i.e., it follows a logarithmic law. 
There has been no analytical results for the degree distribution for this case as, we 
believe, it is extremely difficult if not impossible. 

We now construct the Markov chain for the degree sequence {Ki{t),t = i,i + 1, ...}. 
The state space is = {rui, rm + 1, ...}, where mj = m[lni]. At time t, the probability 
that an existing node i will connect with the new node is given by 



mlnt- 



ki In t 



J2jkj 2t{\nt-l)' 
Hence, the one-step transition probabilities are 



Pkj{t + 1) = P{Ki{t + 1) = i I Ki{t) = k}=< 



1 



2i(lnt-l) ' 



fc In t 
2t(lnt-l) ' 

0, 



(31) 



j = k 

j = k + l 
otherwise 



(32) 



for k = rrii, rrii + t — i, and 



Pkj{t + 1) 



j ^ k 



fov k > rrii + t — i. The transition probability matrix is 



■t rrij In t 

^ 2t(lnt-l) 



rrij In t 
2t(lnt-l) 



1 - 



(mj + l) Int 
2t(lnt-l) 



{rni+l)\nt 
2t(lnt-l) 



(33) 



1 - 



{mi+t—i)\nt {mi-\-t—i)\nt 



2t(lnt-l) 



2t(lnt-l) 



(34) 

{oi t = i,i + 1, ... . 

We note that the structure of the transition matrices here is similar to that of (12). 
The difference is that now m,- is not a constant, in general, but a step function of i, as 
shown in Table 4. 



Table 4: Intervals of rrii keep constant 



time i 


21 


55 


149 


404 


1097 


2981 


8104 


22027 


59875 


162755 


[Ini] 


3 


4 


5 


6 


7 


8 


9 


10 


11 


12 
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Therefore, relations (16), (17) and (18) hold for each interval, e.g., the interval 
(404, 1096). Thus we obtain the following important result 

p(21,t) _ p(21,54) p(55,148) p(59875,162754) p(162755,t) for\ 

— -I- • • • -I- -t*t+i -I- -f t+1 • K'^'^J 

Similarly, the initial probability distribution is fi{i) = (1,0,0, ...) = ei for any i. Thus 
the same algorithm based on (20) can be used to compute the degree distribution P{k, t) 
for this network. 

Prom the computation results, we plot the log — log curves for P(fc, t) for some 
different to and t as shown in Figure 3. We also list some mimerical results in Table 
5. From the figure and the table, it is clear that this network self-organizes into a non- 
stationary scale-free network, with the degree exponent 7 » 3.1 and a positive, though 
very small, non-stationary exponent z. 

Table 5: Numerical results of the logarithmic function case 



parameter m 


time t 


exponent 7 


coefHcient c 


1 


150000 


3.169873 


542.9149 


3 


100000 


3.117526 


1539.876 


3 


150000 


3.081926 


1722.288 


3 


200000 


3.050253 


1952.588 


5 


150000 


3.029171 


2823.081 



4 Conclusions and discussions 

In summary, we introduce a Markov chain-based new method to calculate degree 
distributions of scale- free networks numerically. Comparing with the existing analytical 
methods, this method is more flexible. It offers the asymptotic property of the degree 
distribution for the more complicated models. Using only the transition probability 
matrix Pi{t+1), we can compute the degree distribution P{k). Since the complexity of 
our algorithm is 0(t^), its advantage over the simulation method is also quite obvious: 
it is fast and, for problems that it can handle, it is more reliable and provides better 
understanding of the network behavior. 

The Msc of Markov chain to model the degree evolution is quite novel and opens 
the door for the applications of methodologies and results from a very mature field to 
the exciting new field of scale-free networks. For instance, we may consider to compute 
the joint degree distribution of a node pair by using Markov chains. Furthermore, 
the fact that the evolution of a complex network can be modeled by Markov chains 
may indicate an important direction for us to investigate the underlying mechanisms of 
growth networks, since we have accumiilated extensive understanding of the structural 
properties of Markov chains as we use them to study many natural phenomena. 
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Figure 2: The degree distribution of the power function case 



In Figure El the five lines from left to right correspond to three cases: (1) m = 
1, t = 150,000; (2) TO = 3, t = 100,000, 150,000, 200,000; (3) m = 5, t= 150,000. 
In (2), the three lines are separated, demonstrating the non-stationarity of the degree 
distribution. Again, we can see that the degree exponents are essentially independent 
of TO as the lines are parallel to each other. 




Figure 3: The degree distribution of the logarithmic function case 



In Figure the five lines from left to right correspond to three cases: (1) to = 
1, t = 150, 000; (2) m = 3, t= 100, 000, 150, 000, 200, 000; (3) m = 5, t = 150, 000. In 
(2), the three lines are very close to each other but not entirely overlapping, showing 
that while the degree distribution is not stationary, the non-stationary exponent is very 
small. Again, we can see that the degree exponents are essentially independent of to as 
the lines are parallel to each other. 
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